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Abstract 

We argue that the classical evolution of small x modes in the collision of two ultrarela¬ 
tivistic nuclei is described on a transverse lattice by the Kogut-Susskind Hamiltonian 
in 2-|-1-dimensions coupled to an adjoint scalar field. The initial conditions for the 
evolution are provided by the non-Abelian Weizsacker-Williams fields which consti¬ 
tute the classical parton distributions in each of the nuclei. We outline how lattice 
techniques developed for real time simulations of field theories in thermal equilibrium 
can be used to study non-perturbatively, thermalization and classical gluon radiation 
in ultrarelativistic nuclear collisions. 


1 Introduction 


It is of considerable theoretical and experimental interest to understand the collisions 
of nuclei at ultrarelativistic energies and the putative evolution of the hot and dense 
matter created in these collisions into a thermalized, deconfined state of matter called 
a quark gluon plasma. The theoretical challenge is to understand the dynamics of 
the formation of this matter from QCD, and its properties, while the experimental 
challenge is to detect evidence that such a plasma was indeed formed Q . 

The space-time evolution of the nuclei after the collision and the magnitudes and 
relevance of various proposed signatures of this hot and dense matter depend sensi¬ 
tively on the initial conditions for the evolution, namely, the parton distributions in 
each of the nuclei prior to the collision. In the conventional perturbative QCD ap¬ 
proach to the problem, observables from the collision may be computed by convolving 
the parton distributions of each nucleus, determined from deep inelastic scattering 
experiments, with the elementary parton-parton scattering cross sections. The cross 
sections thereby obtained are often incorporated either in a multiple scattering for¬ 
malism [H or in a classical cascade approach to obtain the space-time evolution |Q. 
While the above approach provides a reasonable description of large transverse mo¬ 
mentum processes at large x, it is not sufficient as we go to small x or alternatively, 
towards central rapidities Q. This is because at small x partons in one nucleus may 
“see” more than one parton in the direction of the incoming nucleus resulting in a 
breakdown of the above described convolution of distributions. What is needed there¬ 
fore to describe the collision of the “wee” nuclear partons is, roughly put, products 
of amplitudes as opposed to products of probabilities. 

A model describing the small x modes in large nuclei was formulated by McLerran 
and Venugopalan |^. The model contains one dimensionful parameter, x(y,Q^), 
which is the total color charge squared per unit area integrated from the rapidity y of 
interest to the beam rapidity. Since it is the only scale in the problem, the coupling 
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constant runs as a function of this scale. One therefore has weak coupling in the 
limits where the color charge x is large; either for ^>>lors^(X). It was argued 
there that the classical background fields in this model are non-Abelian Weizsacker- 
Williams fields. Exact analytical expressions for these fields have been obtained 
recently |^, Further, it has been shown explicitly that x obeys renormalization 
group equations in y and which reduce to the well known DGLAP and BFKL 
equations [^] in the appropriate limits |^, ^|. 

The model was applied to the problem of nuclear collisions by Kovner, McLerran 
and Weigert, who formulated the problem as the collision of Weizsacker-Williams 
fields IQ. The classical background fields after the collision then correspond to solu¬ 
tions of the Yang-Mills equations in the presence of static, random sources of color 
charge on the light cone. The classical background field for the two nuclei after the 
collision was found and perturbative solutions obtained for modes with transverse mo¬ 
menta kt » OiSy/X- After averaging over the Gaussian random light cone sources, 
the energy and number distributions of physical gluons were computed. Further, 
the classical gluon radiation from these perturbative modes was studied by these au¬ 
thors and later in greater detail by several others 1,0 HI- While the perturbative 
approach is very relevant and useful, it is still essential to consider the full non- 
perturbative approach for the following reasons. Firstly, the classical gluon radiation 
computed perturbatively is infrared singular and has to be cut-off at some scale. It 
was argued in Ref. |^, |l^ that a natural scale where the distributions are cut-off 
is given by kt ~ otsy/x- However, since quantitative differences can be large, it is 
important to perform a full calculation. Secondly, the non-perturbative approach is 
crucial to study the possible thermalization of the system and the relevant time scales 
for thermalization. This in turn has several ramifications for computations of various 
signatures of the quark gluon plasma. For instance, if thermalization does occur, then 


as proposed by Bjorken [1T| hydrodynamic evolution of the system is reasonable. In 
that event, our approach would provide the initial temperature and velocity profiles 
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necessary for such an evolution (see |17] and references therein). 

In this paper, we outline how one may perform real-time simulations of the full, 
non-perturbative evolution of classical non-Abelian Weizsacker-Williams fields. Such 
a simulation is possible since the fields are classical. Similar real time simulations 
of classical fields have been performed in the context of sphaleron induced baryon 


number violation \\2\ and chirality violating transitions in hot gauge theories 
In brief, the idea is as follows. We write down the lattice Hamiltonian which describes 
the evolution of these classical gauge fields. It turns out to be the Kogut-Susskind 
Hamiltonian in 2+1-dimensions coupled to an adjoint scalar field. The lattice equa¬ 
tions of motion for the fields are thereby determined straightforwardly. The initial 
conditions for the evolution are provided by the Weizsacker-Williams fields for the 
nuclei before the collisions. Interestingly, the dependence on the static light cone 
sources does not enter through the Hamiltonian but instead from the initial condi¬ 
tions. Also, to reiterate, our results have to be averaged over by the above mentioned 
Gaussian measure. 

The paper is organized as follows. In the following section we briefly discuss 
the problem of initial conditions as formulated by Kovner, McLerran and Weigert 
and their perturbative solution. In section 3, we derive the expression for the lattice 
Hamiltonian. The lattice equations of motion and initial conditions are discussed in 
section 4. We conclude with a discussion of observables that can be computed on the 
lattice and comment on checks that can be performed on such computations. 


2 The classical background field of two nuclei on the 
light cone 

In the work of McLerran and Venugopalan Q, the classical gluon field at small x 
for a nucleus in the infinite frame is obtained by solving the Yang-Mills equations 
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in the presence of a static source of color charge p°‘{rt,rj) on the light cone. Exact 
solutions for the classical field as functions of p°'{rt, rj) were found by Jalilian- Marian 
et al. @ and independently by Kovchegov Distribution functions are computed 
by averaging over products of the classical fields over a Gaussian measure in p with 
the variance Here is the color charge squared per unit area per unit 

rapidity resolved at a scale by an external probe. It is related to x by the 
expression x{v, Q^) = Q^)- 

The above picture of gluon fields in a nucleus at small x was extended to describe 
nuclear collisions by Kovner, McLerran and Weigert |^. We shall review and discuss 
their paper below. The classical background field of two nuclei is described by the 
Yang-Mills equations in the presence of two light cone sources-one on each light cone. 
We have then 


= r , (1) 

where 

= S^+gpl{rt)6{x-) + 6^-gpl{n)6{x+) . (2) 

Gluon distributions are simply related to the Fourier transform A^{kt) of the solution 
to the above equation by < Af{kt)A^{kt) >p. The averaging over the classical charge 
distributions is defined by 

{0)p = j dp+dp-0{p+,p-) 

X (3) 

The averaging over the color charge distributions is performed independently for each 
nucleus with equal Gaussian weight p?. 

The observant reader will notice that we have omitted the rapidity dependence 
of the the charge distributions in the equations immediately above. We will justify 
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this omission in our discussion of the lattice Hamiltonian. We note that the rapidity 
dependence of the charge distribution is also absent in Ref. |^]. 

Before the nuclei collide (t < 0), a solution of the equations of motion is 


= 0 , 

= 6{x~)9{-x^)a\{rt) + 9{x'^)6{-x~)a2{rt), (4) 


where @ 


a 


1,2 


in) = 


1 .dv'-^P±W,rt)'' 


Pe 


V* Pe 


-igl 




, dr]'^ p±{rj' ,rt) 


^9 


( 5 ) 


Above, T] = r]proj — log(x~/x~^gj) is the rapidity of the nucleus moving along the 
positive light cone with the gluon field a\ and rj = —r]proj + log{Xp^gj /x^) is the 
rapidity of the nucleus moving along the negative light cone with the gluon field a^- 
It is expected that at central rapidities (or x << 1) the source density varies slowly 
as a function of rapidity and a* = a^irt)- The above expression suggests that for 
t < 0 the solution is simply the sum of two disconnected pure gauges. 

For t > 0 the solution is no longer pure gauge. Working in the Schwinger gauge 


x+A“+x“A+ = 0, (6) 

the authors of Ref. Q found that with the ansatz 

= ±x^a{T,rt), 

= a\{T,rt), (7) 

where r = \/2x+x“, Eq. ^ could be written in the simpler form 
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] = 0 


rdrT^drO + [Di, D*, a 


;[A, 5^ai] + igT[a, dro] 


-drTdra\ — Z)*, a ] — [D-^ , F-^*] 


0 , 

0 . 


(8) 


The initial conditions for the fields a(r, rt) and given in terms of the fields for 

each of the nuclei at t < 0. We have 


ai|r=o = a\+a2^ 

a|.=o = |[ai,«|]. (9) 


Further, since the equations are very singular at r = 0, the only condition on the 
derivatives of the fields that would lead to regular solutions are c 1 t-q:|t-= 0 ) 9T-Q:^|r=o = 0. 

In Ref. Q], solutions were found in the perturbative limit by expanding the initial 
conditions and the fields in powers p or equivalently, in powers of asp/kt- Performing 
a gauge transformation of the above fields (such that the new fields satisify the 
Coulomb gauge condition), at late times r >> asp, the new fields can be expanded as 
Fourier series with coefficients 2 (^ 4 ) their complex conjugates. The multiplicity 
distribution of classical gluons can then be written as 


dN 

dyd?kt 


1 

(27r)^ 




( 10 ) 


Detailed expressions for classical gluon radiation in the perturbative limit were ob¬ 
tained in Refs. |^, |^. However, all these expressions are infrared singular and 

have to be regulated by an infrared cutoff. One advantage of solving the Yang-Mills 
equations to all orders in asp/h is that it will likely provide a self-consistent, in¬ 
frared safe result for the multiplicity of classical gluon radiation in ultrarelativistic 
nuclear collisions. 
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3 Derivation of lattice Hamiltonian 


While the Yang-Mills equations discussed above can be solved perturbatively in the 
limit askt « it is unlikely that a simple analytical solution exists for Eq. |I| in 
general. The classical solutions have to be determined numerically for t > 0. The 
straightforward procedure would be to discretize Eq. || but it will be more convenient 
for our purposes to construct the lattice Hamiltonian and obtain the lattice equations 
of motion from Hamilton’s equations. 

We start from the QCD action (without dynamical quarks) which is given by 

Sqcd = J , (H) 

where g = det\g^u\- In the forward light cone {t > 0) it is convenient to work with the 
T,r],Ft co-ordinates where r = V‘2x^x~ is the proper time, g = \ \og{x'^/x~) is the 
space-time rapidity and = (x, y) are the two transverse Euclidean co-ordinates. 
In these co-ordinates, the metric is diagonal with g'^'^ = —g^^ = —g^y = 1 and 
= -1/t2. 

After a little algebra, the Hamiltonian can be written as [^] 

H = J dgdrtT | + ^FxyFxy + fAr, + . (12) 

Here we have adopted the gauge condition of Eq. p, which is equivalent to requiring 
A'’" = 0. Also, = ^drAri and = xdrAr are the conjugate momenta. 

Consider the field strength Frfr in the above Hamiltonian. If we assume approx¬ 
imate boost invariance, or 

Ar{T,g,rt) Pi Arir.rt)] Arj{T,g,rt) Pi (13) 

we obtain 

= -Dr^^ , (14) 
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where = dr — igAr is the covariant derivative. Further, if we express in terms 
of the defined in Eq. ^ we obtain the result that = 0 for r > 0. Finally, since 

$ = T^a(r,n); A = a^(r,rl), (15) 


we can perform the integration over the space-time rapidity to re-write the Hamil¬ 
tonian in Eq. as 


H = J dftTT] < - 


1 fda\ 


dr 


^4 xy w 2^2 




/ 


N 2') 

[Dg 





^ J 


(16) 


Here the index f3 = {T,ft). We have thus succeeded in expressing the Hamiltonian 


in Eq. as the Yang-Mills Hamiltonian in 2-|-1-dimensions coupled to an adjoint 
scalar. The discrete version of the above Hamiltonian is well known and is the Kogut- 


Susskind Hamiltonian |]T| in 2-|-1-dimensions coupled to an adjoint scalar field. We 
shall now restrict the discussion to an SU(2) Yang-Mills theory, in order to keep no¬ 
tation simple. Generalization to an arbitrary SU(N) gauge group is straightforward. 
The Kogut-Susskind analogue of Eq. 1 is 


Hl = — 
2t 


EfEt + r 




□ 


1 - -Trf/o 


j,fi j 

where Ei are generators of right covariant derivatives on the group and Uj^fi is a 
component of the usual SU(2) matrices corresponding to a link from the site j in the 
direction h. The first two terms correspond to the contributions to the Hamiltonian 
from the chromoelectric and chromomagnetic field strengths respectively. In the last 
equation is the adjoint scalar field with its conjugate momentum p = p°‘(T°‘. 

Finally, we should comment on a key assumption in the above derivation, namely, 
the boost invariance of the fields. This invariance results in Eq. |l^ thereby allowing us 
to restrict ourselves to a transverse lattice alone. To clarify the issue we are compelled 
to make a few historical remarks. As we mentioned earlier, the authors of Ref. Q 
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found a solution which was explicitly boost invariant. However, this result was a 
consequence of the original assumption of McLerran and Venugopalan that the color 


charge density factorizes, p°'{rt,ri) —> p^{rt)6{x~). It was noticed in Ref. |^], this 
factorized form for the charge density results in infrared singular correlation functions 
which diverge as the square of the lattice size. This problem was resolved in Ref. Q 
where the authors realized that a rapidity dependent charge density p'^{rt,r]) would 
give infrared safe solutions. This might be interpreted as implying that the boost 
invariance assumption of Ref. should be given up as well. 

Fortunately, this is not necessary. In principle, the rapidity dependence of the 
color charge density can be arbitrarily weak since that is sufficient to obtain infrared 
safe correlation functions. In Ref. [21|, an explicit model was constructed for the 
color charge distribution in the fragmentation region. It was shown there that for 
V < Vproj the color charge distribution had a very weak dependence on rj. We should 
note too that it was shown recently by Gyulassy and McLerran Q that the initial 
conditions in Eq. 0 are unaffected by the smearing in rapidity. 


4 Lattice equations of motion and initial conditions 


Lattice equations of motion follow directly from Hi oi Eq. H- Eor any dynamical 
variable v with no explicit time dependence v = where v is the derivative 

with respect to r, and {} denote Poisson brackets. We take Ei, Ui, pj, and Aj as 
independent dynamical variables, whose only nonvanishing Poisson brackets are 


K, ^>5} = 6,j5ab] {Ef, Um} = {Ef, E^} = 25i^eabcEf 


(no summing of repeated indices). The equations of motion are consistent with a set 
of local constraints (Gauss’ laws). These are 


h 



a 


ITT 



- - 


2eabcp]A^j = 0 . 


(18) 
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The initial conditions for the helds and momenta in Eq. ^ can be discretized as 
follows. To every lattice site j we assign two SU(2) matrices, Vij and V 2 ,j- Each of 
these two defines a lattice gauge field conhguration (corresponding to the two pure 
gauges discussed in section 2) whose link variables are 

where g = 1, 2 labels the two nuclei. An initial condition for the link matrices, having 
the correct formal continuum limit, is then 


^j',ri|r=0 


-1 




1/2 ’ 


(19) 


where I is the identity matrix. Introducing Uq through exp(—faag) = one can 
easily verify that for smooth initial held conhgurations 

1 

1 - -Un —B^B^ 

2 2 

in the limit of vanishing lattice spacing a. Here B'^ is the i^th color component of the 
SU(2) magnetic held corresponding to the gauge potential ai + a 2 - Thus the initial 
plaquette term of Eq. 0 has the correct formal continuum limit. 

As we have already discussed, regularity requires the vanishing of the initial 


transverse color electric helds. If follows from Eq. 15 that the adjoint scalar A must 
also vanish initially. On the other hand, values of p, the conjugate momentum of A, 
need not be zero at r = 0. In the continuum 


P^\t=o = 2a. 


This initial condition can be written on the lattice as follows: 

Pj\r=0 = + P]-h,hM‘j-h,h^ ■ 

h 
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Note that the initial conditions as described satisfy the Gauss constraints [|. 

Finally, we comment on numerical integration of the lattice equations of motion. 
A method of choice for energy-conserving Hamiltonian equations is the leapfrog al¬ 
gorithm, wherein, the fields are updated at every step of the integration while their 
conjugate momenta are kept fixed and vice versa. In our case, the Hamiltonian de¬ 
pends explicitly on time and the energy is not conserved. Nevertheless, only a minor 
modification of the algorithm is necessary in order to maintain the same time step ac¬ 
curacy as in the energy-conserving case. As explained in the leapfrog algorithm 
has the useful property of respecting Gauss’ laws exactly (regardless of the time step) 
for any Hamiltonian which is a sum of potential and kinetic terms. 

5 Outlook 

We have outlined above a procedure to solve for the gauge field configurations pro¬ 
duced in the collision of the Weizsacker-Williams fields of two nuclei. These can 
then be used to compute a large variety of observables as a function of proper time 
r. An observable that can be computed directly in our Hamiltonian approach is 
the energy density (and correlations in the energy density). Also, by looking at the 
Fourier decomposition of the co-ordinate space correlators of the electric fields, one 
can determine the energy and number distributions of the modes as well as the energy 
dispersion relation for the soft modes. These can be directly related to the multiplic¬ 
ity and transverse energy of mini-jets (for a review see Ref. |23||). Some care however 
must be exercised in ensuring that the residual gauge freedom of the fields is fixed 
properly. 

The field configurations generated immediately after the collision are completely 
^The lattice initial conditions are not completely specified by requiring that they reduce to the 
corresponding conditions in the continuum upon taking the formal continuum limit. This arbitrari¬ 
ness can be removed by a consistent derivation of the initial conditions entirely within the lattice 
theory. Such derivation will be presented elsewhere p^. 
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out of equilibrium. It is an interesting problem to study if and how the system 
approaches equilibrium. This can be done by studying whether the multiplicity 
saturates as a function of time. Alternatively, one can follow the example of the 
Duke group and study the behaviour of the Lyapunov exponents as a function of 
time 1^, If the system does indeed thermalize, energy density/temperature 

and velocity profiles can be extracted for use as initial conditions in hydrodynamic 
simulations. 

A useful check of the results of our proposed simulation (besides the usual tech¬ 
nical ones) is to reproduce the results of Ref. I, I 0] in the “Abelian limit” of large 
transverse momenta, kt» asfJ-- The sensitivity of various observables to the lattice 
spacing and lattice size must also be carefully studied. 

There are several open questions which have still to be resolved. Primarily, our 
simulation is completely classicaRhow accurate is the classical description? It would 
be interesting to study whether the Wilson renormalization group evolution for the 
scale x(ry, can be implemented on the lattice. Another issue is our assumption of 
boost invariance-relaxing this assumption however is not conceptually difficult but 
may be numerically time consuming. 

A more detailed discussion of the above and detailed numerical simulations will 


be presented at a later date [ 261. 
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